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Abstract 

A simple closed-form analytic expression for the probability of two-flavour neutrino 
oscillations in a matter with an arbitrary density profile is derived. Our formula is 
based on a perturbative expansion and allows an easy calculation of higher order cor- 
rections. The expansion parameter is small when the density changes relatively slowly 
along the neutrino path and/or neutrino energy is not very close to the Mikheyev- 
Smirnov-Wolfenstein (MSW) resonance energy. Our approximation is not equivalent 
to the adiabatic approximation and actually goes beyond it. We demonstrate the va- 
lidity of our results using a few model density profiles, including the PREM density 
profile of the Earth. It is shown that by combining the results obtained from the 
expansions valid below and above the MSW resonance one can obtain a very good 
description of neutrino oscillations in matter in the entire energy range, including 
the resonance region. 
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1 Introduction 



In most neutrino oscillation experiments neutrinos propagate substantial distances in mat- 
ter before reaching a detector, and therefore an accurate description of neutrino oscillations 
in matter [1, 2] is an important ingredient of the analyses of the data. For a matter of an 
arbitrary density profile the neutrino evolution equation admits no closed-form solution, 
and one usually has to resort to numerical methods. While numerical integration of the 
evolution equation usually poses no problem, it is still highly desirable to have approximate 
analytic solutions, which may provide a significant insight into the physics of neutrino os- 
cillations in matter, clarify the dependence of the oscillation probabilities on the neutrino 
parameters and in many cases help save the CPU time. To this end, a number of analytic 
solutions of the neutrino evolution equation in matter, based on various approximations, 
has been developed (for recent studies, see e.g. [3, 4, 5, 6, 7, 8, 9, 10]). 

In this paper we derive a simple analytic expression for the two-flavour oscillation prob- 
ability valid for an arbitrary matter density profile. We employ a perturbative approach 
based on the expansion in a parameter which is small when the density changes relatively 
slowly along the neutrino path and/or neutrino energy is not very close to the Mikheyev- 
Smirnov-Wolfenstein (MSW) [1, 2] resonance energy. Our approximation is not equivalent 
to the adiabatic approximation and actually goes beyond it. We demonstrate the validity 
of our results using a few model density profiles, including the important PREM profile 
[11], which gives a realistic description of matter density distribution inside the Earth. We 
also show that, by combining the results obtained for the energies below and above the 
MSW resonance ones, one can obtain an excellent description of neutrino oscillations in 
matter in the entire energy range. The simple form of our result and the wide range of its 
applicability are the two main advantages of this approach. 

An approach similar to ours has been employed in [12, 7]. Unlike in those publications, 
in the present work we do not confine ourselves to the leading approximation, but also 
calculate the first and second order corrections and show that this improves the accuracy 
of the approximation drastically. 

The paper is organized as follows. In Sec. 2 we present the formalism used to derive our 
analytic solution. In Sec. 3 we apply this method to the case of a parabolic and a power 
law matter potentials. In Sec. 4 we present the results obtained in the case of the realistic 
PREM Earth's density profile. We discuss our results and conclude in Sec. 5. 

2 The formalism 

In a number of important cases the full three-flavour neutrino oscillations can to a very 
good accuracy be reduced to effective two-flavour ones. These include v e <-> Vy,(v T ) os- 
cillations either in the limit of vanishingly small 1-3 mixing, when the oscillations are 
essentially driven by the "solar" parameters Am^ and #12, or at sufficiently high energies 
(E > 1 GeV for oscillations in the Earth), when the 1-2 mixing in matter is strongly sup- 
pressed; in that case the oscillation probabilities are essentially independent of the "solar" 
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(1) 



parameters and are governed by Am| x and #13. For definiteness, in our numerical examples 
we will concentrate on the second case, though our general discussion will be valid in both 
situations. 

Two-flavour oscillations of neutrinos in matter are described by the Schrodinger-like 
evolution equation [1, 2] 

t) = ( B Aj(l 

where the overdot denotes the differentiation with respect to the coordinate, and £ and 
rj are respectively the probability amplitudes to find v e and u a , the latter being a linear 
combination of and v T . In the limit when the 1-3 mixing vanishes, #13 — > 0, one has v a = 
cos#23^V — sm @23V T , whereas in the situations when the solar parameters play practically no 
role (e.g. for oscillations of high-energy neutrinos in the Earth), v a = sin #23^ + cos023^r- 
The quantities A and B in Eq. (1) are 

B = 5 sin 29 , 
A(x) = 8 cos26 -V(x)/2. (2) 

Here the function A(x) depends on the electron number density N e (x) through the Wolfen- 
stein potential V(x) defined as 

V(x) = V2G F N e (x) = 7.54 x 10" 14 Y e (x) p(x)(g/cm 3 ) eV, 

where Gp is the Fermi constant, p(x) is the mass density of matter and Y e (x) is the number 
of electrons per nucleon. The parameter 5 is defined as 5 = Am 2 /4£, and 8 is the relevant 
mixing angle in vacuum. In the limit #13 — > one has Am 2 = Am|i, #0 — $12 , and the 
v e ^ y^iy-i) oscillation probabilities are given by 

P{v e -> x) = Piy^ -> V e ] x) = cos 2 9 23 P 2 (x) , (3) 

P(u e -> v T yx) = P[y T -> i/ e ;x) = sin 2 6 23 P2 (x) . (4) 
Here P2{x) is the effective two-flavour oscillation probability: 

P 2 {x) = P(u e ^u a ;x) = |t/(x)| 2 (5) 

(we assume the initial conditions £(0) = 1, 7/(0) = 0). For oscillations of high-energy 
neutrinos in the Earth one has Am 2 = Am 2 1; #0 = $13> and the ^ e <-» v^{v T ) oscillation 
probabilities are 

^(^ e -> f M ; a?) = P^y ji -> Ve\ x) = sin 2 6*23 , (6) 

P(u e -> i/ r ; x) = P(r/ r -> i^ e ; sc) = cos 2 6» 23 ^2(2) , (7) 
where, as before, P2{x) is given by Eq. (5). 
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Differentiating Eq. (f), one can find decoupled second order differential equations for 
and rj(x) [13, 14]. The equation for the transition amplitude r/(x) reads 

f] + (u 2 + iA)r) = , (8) 

where we have defined the function uj{x) as 

uo\x) = A 2 (x) + B 2 . (9) 

Note that the instantaneous eigenvalues of the effective Hamiltonian in Eq. (1) are ±u(x). 
The equation for differs from Eq. (8) by the sign of the A term. 

It will be convenient for our purposes to rewrite Eq. (8) in the following form: 

fj + (uj 2 - iu)ij = (-iA)77, (10) 

where we have introduced the notation 

A = A + u. (11) 

Eq. (10) cannot in general be solved exactly, but, as we shall see, it admits a simple 
perturbative solution. To show that, let us first notice that, for energies (or densities) 
above the MSW resonance one, the quantity A on the right hand side of Eq. (10) is small. 
Indeed, from Eqs. (9) and (2) it follows that for V/2 — cos 20 o 5 ^> sin2# < i > (i-e. for 
—A 3> B) one has <jj ~ —A, so that A ~ 0. The smallness of the parameter A allows one 
to solve Eq. (10) perturbatively, order by order. Expanding in powers of A, we find the 
equation for the nth order transition amplitude r) n (with n > 0): 

fj n + (lj 2 - iu)r} n = (-zA)r/ n _i . (12) 

The zero order transition amplitude r] satisfies the equation with the vanishing right hand 
side: 

rjo + (u 2 - iu)r] = . (13) 

Its solution for an arbitrary functional dependence of u(x) on the coordinate can be readily 
found by considering the quantity X = rj — iur] , which, as follows from (13), satisfies 
the first-order equation X + iuoXq = 0. Taking into account that the initial conditions 
£(0) = 1, 77(0) = also imply, through Eq. (1), r)(0) = — iB, one finds 



where 



Be* ] f dx 1 e- 2i ^ Xl \ (14) 
Jo 



x 



uj{x')dx' . (15) 
This yields the zero-order solution for the two-flavour transition probability P2{x) [12, 7]: 



[P 2 {x)] = \vo(x)\ 2 = B 



2 

dxi e -2ut>( X1 ) 



(16) 
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Assuming that the amplitude r] n _i(x) on the right hand side of Eq. (12) is known, one 
can solve it for rj n . To this end, we introduce the quantity 

X n = Ti n - iwq n , (17) 

in terms of which Eq. (12) can be rewritten as 

X n + icoX n = (-iA)^-! . (18) 

This can now be solved by the standard methods. First, we find the general solution of 
the homogeneous equation 

X n + iiuX n = 0, (19) 

which gives 

X n (x) = F e~ l ^ x) (20) 

with F an integration constant. Next, the solution of the inhomogenous equation (18) is 
found by allowing F to depend on the coordinate x and substituting Eq. (20) back into 
Eq. (18). Taking into account the initial condition F(0) = fj(0) — iu(0)r)(0) = —iB, one 
finds 

F(x) = J dx x e^ {xi) (-iA (zi)) Vn-i (zi) - iB . (21) 

The solution for X n is now given by Eq. (20) with F replaced by F(x) from Eq. (21). Once 
X n is known, it is straightforward to solve Eq. (17) for r) n . This yields 

T] n (x) = e***) f d Xl e~ 2i ^ r dx 2 e*™ (-iA (x 2 )) r^i (x 2 ) + r) {x) , (22) 



where we have used Eq. (14). The corresponding nth order effective two-flavour oscillation 
probability is then found as [P 2 (^)]n = \Vn(x)\ 2 . 

Eq. (22) represents the main result of our paper. It gives an analytic expression for the 
oscillation amplitude in the nth order in perturbation theory in terms of the lower-order 
solutions Tjn-i and Tjo- For our numerical illustrations we will consider the solutions with 
n = 0, 1 and 2. 

Eq. (22) has been derived under the assumption that A is a small parameter. As we 
pointed out before, this is true for energies above the MSW resonance one. This means 
that the perturbative approach considered above should, in general, fail for energies below 
the MSW resonance one. However, a simple modification of the above procedure leads to 
a description of neutrino oscillations valid below the MSW resonance. In order to show 
this, let us, instead of casting Eq. (8) in the form (10), rewrite it as 

f) + (u 2 + iuj)t} = (-iA)T) , (23) 

where A is now defined as 

A = A - Co . (24) 
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For small vacuum mixing angles, this is a small parameter below the MSW resonance, 
since in that case A>B and so u ~ A. Therefore, we can proceed with the perturbative 
approach, as before. Comparing Eqs. (23) and (24) with Eqs. (10) and (11) respectively, we 
see that the two pairs of equations differ only by the sign of u(x). Therefore the solution 
of Eq. (23) can be obtained from Eq. (22) by simply replacing uj(x) by —uj(x). This will 
also change the values of the oscillation probabilities obtained in all orders in perturbation 
theory except for the zero-order probability which, as can be seen from (16), is invariant 
with respect to the flip of the sign of u(x). As we shall see, by combining the results valid 
above and below the MSW resonance one can obtain a very good description of neutrino 
oscillations in matter in the entire energy range. 

Let us now discuss the expansion parameter of our perturbative approach. We have 
found that the corrections to the zero order amplitude 770 are proportional to A = A ± u, 
where the upper and lower signs refer to the energies above and below the MSW resonance, 
respectively. These quantities can be expressed through the mixing angle in matter Q Tx •■ 



V 

A = A±lu = — — [1 ± cos20 m ] . (25) 

Far above the MSW resonance one has cos2# m ~ —1, whereas far below the resonance 
cos 29 m ~ cos 26*o, which is close to 1 in the case of small vacuum mixing. This demonstrates 
the smallness of A in its corresponding domains of validity. At the MSW resonance one 
has cos26* m = 0, and A is only small if V is. 

An examination of Eq. (22) shows that the expansion parameter of our perturbative 
approach is actually ~ |A|/kj 2 (see Eq. (15)). In various energy domains we have 
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if (c2<5 — V/2) ^> S2S (below the resonance) 
if |c2<5 — V/2 1 s 2 5 (near the resonance) 
if (V/2 — c 2 5) ^> s 2 S (above the resonance) 

(26) 

where we have used the shorthand notation c 2 = cos 28q, s 2 = sin26*o. From Eq. (26) it is 
easy to see that outside the MSW resonance region the expansion parameter approximately 
satisfies 

j^- ~ sin 2 20 m ^ = sin 29 m 7^ , (27) 

where 7msw = 4:U 3 / (\V\B) = 4u 2 /(\V\ sin2# m ) is the MSW adiabaticity parameter. Thus, 
for small mixing in matter (sin2# m <C 1) our approximation is better than the adiabatic 
one. Close to the resonance the two approaches have comparable accuracy. 



L Note that sin2# m = B/u, cos2^ TO = A/uj. 
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3 Two examples: parabolic and power law profiles 



As a first study, we apply our formalism to two simple density distributions: a parabolic 
and a power law profile. 

For the parabolic profile, we consider the following density distribution: 



PO) = Po 

with 



- k 



L 2 /4 



(2? 



Po = Pmax = 8 g/cm 3 , k=l--^- = 0.5, (29) 

Pmax 

and we take the baseline to be L = 10000 km. Note that the parabolic density profile 
represents a good approximation for the density distribution felt by neutrinos in the Earth 
when they cross only the Earth's mantle. 

Next, we analyze the case of the following power-law density distribution: 

p(x)=p (^-) 3 (30) 

with 

x = 10 3 km and p = 10 3 g/cm 3 , (31) 

and we consider neutrino propagation over the distance L = 100 km. The profile p oc x~ 3 
represents a realistic description of the density distribution inside supernovae; note, how- 
ever, that neutrino flavour transitions in supernovae are more adequately described by 
different methods (see, e.g., [15]), and so we consider the profile (30) just for illustration. 

The results based on our perturbative analytic approach for the profiles (28) and (30) 
are presented in Fig. 1, where they are compared with the exact ones, obtained by direct 
numerical integration of the neutrino evolution equation (1). The upper panels show the 
oscillation probabilities for the parabolic density profile and the lower ones, for the power- 
law profile (30). The left panels correspond to the expansion valid for energies below 
the MSW resonance ones, whereas the right panels were obtained for the expansion valid 
above the resonance energies. As expected, the zero-order approximation gives a good 
accuracy only outside the MSW resonance region (i.e., outside the intervals E ~ 3 - 6 
GeV for the parabolic profile and E ~ 30 - 50 MeV for the power-law one). 2 The first- 
order perturbative results obtained using the expansion valid below the MSW resonance 
extend slightly the region of good accuracy towards higher energies, closer to the MSW 
resonance, though in general fail for energies above the MSW resonance, whereas the first- 
order results found from the expansion valid above the MSW resonance extend the region 
of good accuracy to lower energies, but in general fail below the MSW resonance. Thus, 
the first-order calculation taken in their respective domains of applicability allow to achieve 



2 Note that, since the profiles (28) and (30) (as well as the PREM profile considered in the next section) 
span a range of matter densities, neutrinos in an interval of energies experience the MSW resonance. 
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a good description of the exact results closer to the resonance energy than the zero-order 
solutions do, i.e. they reduce the energy domain in which the approximation fails. At 
the same time, as can be seen from Fig. 1, the second-order probabilities |t7 2 | 2 practically 
coincide with the corresponding exact results, irrespectively of whether they are obtained 
using the expansion valid below or above the MSW resonance. 

4 Propagation inside the Earth: PREM profile 

Neutrinos coming from various sources can propagate inside the Earth before reaching a 
detector. Examples are atmospheric neutrinos, neutrinos coming from WIMP annihilation 
inside the Earth or the Sun, as well as neutrinos studied in long-baseline accelerator exper- 
iments. We will consider here oscillations of high-energy neutrinos in the Earth, for which 
we take the matter density distribution as described by the PREM profile [11] (Fig. 3). 
Note that the PREM profile is symmetric with respect to the midpoint of the neutrino 
trajectory, and therefore the two-flavour transition amplitude r)(x) obtained as a solution 
of Eq. (1) is pure imaginary due to the time reversal symmetry of the problem [16]. 

In Fig. 2 we present the oscillation probability P2 as a function of neutrino energy 
E for two values of the zenith angle of the neutrino trajectory: cos6 z = —1, when the 
neutrinos propagate the longest distance inside the Earth, traversing it along its diameter, 
and cos #2 = —0.95, when they do not cross the inner core of the Earth. As in Fig. 1, 
we compare the approximate solutions, up to the second order (^2)2 = l^l 2 ), with the 
exact solutions found by direct numerical integration of the neutrino evolution equation. 
In this figure (as well as in Figs. 4-6 below) in the left panels we present the oscillation 
probabilities obtained with the expansion valid below the MSW resonance energy, whereas 
the right panels show the results found from the expansion valid above the MSW resonance. 

It can be seen from Fig. 2 that the zero order probability |?7o| 2 reproduces accurately 
the exact one, \r]\ 2 , only for energies that are outside the resonance region. Indeed, the two 
solutions nearly coincide for E < 2.5 GeV and for E > 7 GeV, but deviate substantially 
between these energies. As can be seen from the figure, the accuracy of the first order solu- 
tions is slightly better in their respective domain of validity: the solutions |?7i| 2 valid below 
the resonance (left panels of the figure) allow an accurate description of the probability 
for slightly higher energies than |?7o| 2 does, allowing to come closer to the MSW resonance 
from below; however, they fail badly (not even being bounded by 1) above the resonance. 
Likewise, the solutions \r]i\ 2 valid above the resonance (right panels) allow to come closer 
to the MSW resonance from above, but fail below the resonance. 

At the same time, the second-order solution \r]2\ 2 gives quite a good approximation to 
the exact probability |?7| 2 for all energies, though the solutions obtained through the ex- 
pansions in their corresponding domains of validity give a better accuracy in these energy 
domains. By combining the second-order solutions valid below and above the MSW reso- 
nance, one can obtain a very good description of the exact oscillation probability practically 
at all energies, including the resonance region. We have also checked that for trajectories 
that do not cross the core of the Earth (cos# z > -0.838), for which the matter density pro- 
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file seen by the neutrinos is relatively smooth, the second order solutions obtained through 
both expansions essentially coincide with the exact one for all energies. 

In Fig. 4 we present the oscillation probability P2, obtained in different orders in per- 
turbation theory, as a function of the distance travelled by neutrinos inside the Earth for 
vertically up-going neutrinos (cos 6^ = —1) and two values of neutrino energy, E = 2.8 GeV 
and 6 GeV. The figure clearly demonstrates how the accuracy improves with increasing 
order in perturbation theory; the second order solutions \r)2\ 2 nearly concide with the exact 
probability \r]\ 2 along the entire neutrino path. 

Fig. 5 illustrates the dependence of the analytic solutions on the zenith angle for two 
values of the neutrino energy, E = 2.5 GeV and 6 GeV. For both energies we show the 
solutions obtained using the expansions valid below and above the MSW resonance. The 
results agree with our expectations: the second order solution based on the expansion 
valid below the resonance reproduces the exact one extremely well for E = 2.5 GeV but 
does not give a good accuracy (especially in the core region) for E = 6 GeV, while the 
situation is opposite in the case of the solution corresponding to the expansion valid above 
the resonance. 

Finally, in Fig. 6 we show the dependence of the accuracy of the analytic solutions on 
the value of the vacuum mixing angle 6q = #13. As one can see by comparing the upper 
panels with the corresponding lower ones, with decreasing value of #13 the accuracy of our 
perturbative expansion improves. This is the consequence of the fact that the expansion 
parameter (27) decreases with decreasing #13. 

5 Discussion and conclusions 

We have developed a perturbative approach for two-flavour neutrino oscillations in matter 
with an arbitrary density profile. The zero-order oscillation amplitude 770 satisfies the 
equation which can be solved analytically for an arbitrary dependence of the matter density 
distribution on the coordinate along the neutrino path; higher order amplitudes are then 
obtained from the lower-order ones through a simple perturabative procedure. We have 
studied the zeroth, first and second order solutions and compared them with each other 
and with the exact oscillation probability obtained by numerical integration of the neutrino 
evolution equation. In all orders except the zeroth one, the expansion scheme depends on 
whether the neutrino energy is above or below the MSW resonance energy, and one has to 
consider these two cases separately. 

While the zero-order result gives a very good accuracy outside the resonance region, 
higher order corrections are necessary to achieve an accurate description of the oscillation 
probability in the vicinity of the MSW resonance. We have demonstrated how these cor- 
rections, when taken in their respective energy domains of validity, improve drastically the 
precision of the approximation. 

For the smooth density profiles that we have studied, we found that the second order 
oscillation probability reproduces the exact one extremely well in the whole interval of 
energies, including the MSW resonance region, independently of whether the expansion 
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scheme valid below or above the resonance was used. The same is also true for the PREM 
density profiles in the case when neutrinos cross only the mantle of the Earth, since the 
density jumps experienced by neutrinos in that case are relatively small. The high accuracy 
of the second order approximation for smooth density profiles is related to the fact that our 
expansion parameter, Eq. (27), is proportional to \V\. This parameter is smaller than the 
expansion parameter of the adiabatic approximation by the factor sin20 m and therefore 
our approach gives a better accuracy than the adiabatic expansion when the mixing in 
matter is small. Note that a different expansion of the same evolution equation (8) was 
employed in [17]. 

For energies above the MSW resonance, our expansion parameter is essentially 

M ~ sin 2 2# m ^. (32) 
or v z 

For the PREM density profile of the Earth, the function |^|/V^ 2 is plotted in the right panel 
of Fig. 3. As can be seen from the figure, in most of the coordinate space the value of 
this function does not exceed 0.25. The spikes corresponding to the density jumps, though 
quite high, are very narrow; they do not destroy our approximation because their contribu- 
tions get suppressed due to the integrations involved in the calculation of the higher-order 
corrections to the oscillation amplitude (see Eq. (22)). Still, these contributions are not 
negligible, especially for neutrinos crossing the Earth's core. As a result, for core-crossing 
neutrinos with energies close to the MSW resonance ones, even the second-order oscillation 
probabilities are only adequate when taken in their respective energy domains of validity. 
By combining the solutions valid below and above the MSW resonance one obtains a very 
accurate description of neutrino oscillations in matter in the entire energy range. 

To conclude, we have derived a simple closed-form analytic expression for the prob- 
ability of two-flavour neutrino oscillations in a matter with an arbitrary density profile. 
Our formula is based on a perturbative expansion and allows an easy calculation of higher 
order corrections. We have applied our formalism to a number of density distributions, 
including the PREM density profile of the Earth, and demonstrated that the second-order 
approximation gives a very good accuracy in the entire energy interval. 
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Figure 1: Oscillation probability P2 versus neutrino energy E in the case of the parabolic 
(upper plots) and power law (lower plots) density profiles. Left panels: probabilities ob- 
tained from the expansion valid below the MSW resonance, right panels: the same for the 
expansion valid above the resonance. We have taken Am 2 = 2.5 10~ 3 eV 2 and Y e = 0.5. 
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Figure 2: Probability P2 versus neutrino energy E for neutrino oscillations in the Earth 
(PREM density profile) for two values of the zenith angle. Left panels: probabilities 
obtained from the expansion valid below the MSW resonance, right panels: the same for 
the expansion valid above the resonance. We have taken Am 2 = 2.5 10 -3 eV 2 . 
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Figure 3: Left panel: matter density distribution inside the Earth as predicted by the 
PREM profile [11]. Right panel: the function —V/V 2 as calculated with the PREM profile 
with density jumps smoothed over the distance of 30 km. 
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Figure 4: Oscillation probability P2 in different orders in perturbation theory versus the 
distance travelled by neutrinos inside the Earth, for cos8 z = -1.0 and for two values of 
neutrino energy (E= 2.8 GeV and 6 GeV). Left panels: probabilities obtained from the 
expansion valid below the MSW resonance, right panels: the same for the expansion valid 
above the resonance. We have taken Am 2 = 2.5 10 -3 eV 2 . 
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Figure 5: Oscillation probability P2 versus the zenith angle 9 Z for neutrinos propagating 
inside the Earth, for neutrino energies E = 2.5 and 6 GeV. Left panels: probabilities 
obtained from the expansion valid below the MSW resonance, right panels: the same for 
the expansion valid above the resonance. We have taken Am 2 = 2.5 10 -3 eV 2 . 
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Figure 6: Oscillation probability P2 versus neutrino energy E in the case of cos #2 = —0.95 
and for two different values of sin 2 29 1^. Left panels: probabilities obtained from the 
expansion valid below the MSW resonance, right panels: the same for the expansion valid 
above the resonance. We have taken Am 2 = 2.5 10 -3 eV 2 . 
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